
function [obj_cqr,gamma0_cqr,gammap_cqr,res_cqr]=CQR(q, p, y, x, toler,gamma0,gammap)
 
 % Returns the Composite Quantile Regression (CQR) estimator for a quantile
 % regression model with constant slope coefficients
 % The CQR estimator is computed using a MM (Majorize-Minimize) algorithm
 % from Hunter Lange (2000) JCGS
 
 %******************* Input Variables **************************************

 % q : quantile level
 % p : dimension of the covariate (including intercept)
 % y : dependent variable
 % x : covariate (must include intercept)
 % toler : tolerance parameter for the precision of the estimation
 % gamma0 : initial value for the intercept 
 % gammap : initial values for the slope coefficients (vector of dimension p-1)
 
 %******************* Output Variables **************************************
 
 % obj_cqr : optimized composite quantile regression objective function
 % gamma0_cqr : estimate of the intercept coefficient
 % gammap_cqr : vector storing the estimates of the slope coefficients
 % res_cqr : residuals of the estimation
 
 
  n=length(y);
  
  %******** MM tolerance and precision variables ****************************

  tn=toler/n;
  e0=-tn/log(tn);
  epsilon=(e0-tn)/(1+log(e0)); %  epsilon is the smoothing parameter.  
  % As suggested in Hunter Lange paper, it
  %  is chosen to roughly satisfy -epsilon(log epsilon)=toler/n using a
  %  Newton-Raphson step above.
  
  %******** Algorithm initialisation *******************
  
  iteration=0;
  change=realmax;
  m=size(q,2);
  for k=1:m
      res(:,k)=y-x(:,1)*gamma0(k)-x(:,2:p)*gammap; % Compute residuals for current gamma for each quantile level k.
      weights(:,k)=1./(epsilon+abs(res(:,k)));
      v(:,k)=1-2*q(:,k)-res(:,k).*weights(:,k);
      lastsurrk(k,:)=-sum(res(:,k).*v(:,k))-sum((1-2*q(:,k)).*res(:,k)); % Last value of surrogate (up
                                             % to a constant).
  end
  lastsurr=sum(lastsurrk)/m; % average of the surrogate over the quantile levels 
   
   %******** Algorithm loop *******************  
  while change > toler 
	iteration = iteration + 1;
	J=x;  % J is the first differential matrix.
    
        for k=1:m
        weightsk=weights(:,k);
       	Hk(:,:,k) = J'*(weightsk(:,ones(p,1)).*J);
        Hp=sum(Hk,3);
        H0k(:,:,k)=Hk(1,:,k);
        Hpk(:,:,k)=[Hk(2:p,1,k),Hp(2:p,2:p)];
        H(:,:,k)=[H0k(:,:,k);Hpk(:,:,k)];
        Dk(:,k) = J' * v(:,k);
        Dp=sum(Dk,2);
        step(:,:,k)=-inv(H(:,:,k))*[Dk(1,k);Dp(2:p,:)];
        end
        step0=step(1,:,:);
        step0=reshape(step0,1,m);
        stepp=step(2:p,:,:);
        stepp=reshape(stepp,p-1,m);
        stepp=sum(stepp,2)/m;
	    gamma0=gamma0+step0;
        gammap=gammap+stepp;
        for k=1:m
       	res(:,k)=y-x(:,1)*gamma0(k)-x(:,2:p)*gammap;
        surrk(k,:) = sum(res(:,k).^2.*weights(:,k)) + sum((4*q(:,k) - 2).*res(:,k));
        end
        surr=sum(surrk)/m;
        
        %  Now do the step-halving procedure to ensure a decrease in the
        %  value of the surrogate function.
        while surr > lastsurr
       		step0=step0/2;
            stepp=stepp/2;
       		gamma0=gamma0-step0;
            gammap=gammap-stepp;
            for k=1:m
       		res(:,k)=y-x(:,1)*gamma0(k)-x(:,2:p)*gammap;
            surrk(k,:) = sum(res(:,k).^2.*weights(:,k)) + sum((4*q(:,k) - 2).*res(:,k));
            end
            surr=sum(surrk)/m;
       	end

       	    change=lastsurr-surr;
            for k=1:m
	        weights(:,k) = 1./(epsilon+abs(res(:,k)));
            v(:,k)=1-2*q(:,k)-res(:,k).*weights(:,k);
	        lastsurrk(k,:)=-sum(res(:,k).*v(:,k))-sum((1-2*q(:,k)).*res(:,k));
            end
            lastsurr=sum(lastsurrk)/m; 
  end
  gamma0_cqr=gamma0;
  gammap_cqr=gammap;
  res_cqr=res;
  for k=1:m
      resk=res(:,k);
      obj(k,:)=sum(q(:,k).*resk)-sum(resk(resk<0));
  end
  obj_cqr=obj;
  